####################################################
# JL and RT. "Testing models of legislative decision-making with measurement error: the robust predictive power of bargaining models over procedural models" #
# Replication Log File #


##############
## FIGURE 2 ##
actual <- rnorm(500, 50, 10)
mod1 <- rnorm(500, 50, 10)

abs.error1 <- abs(mod1-actual)
mean(abs.error1) #11.07

mod2 <- actual-20
abs.error2 <- abs(mod2-actual)
mean(abs.error2) # 20

res1 <- lm(actual~mod1+mod2-1)
summary(res1)


#jpeg("Figures/Figure-2.jpg", width=600, height=800)
nf <- layout(matrix(c(0,1,1,1,1,1,0,0,2,2,2,2,2,0,0,0,3,3,3,0,0), 3,7, byrow=T))

hist(actual, main="",xlim=c(0,100), ylim=c(0,60), xlab="", col="gray40", breaks=35, cex.axis=2.2, ylab="")
mtext("Frequency", 2, line=3.5, cex=1.5)
legend("topright", legend=c("Actual"), fill=c("gray40"), cex=2.5)

hist(mod1, main="",xlim=c(0,100), ylim=c(0,60), xlab="", breaks=35, col="white", cex.axis=2.2, ylab="")
mtext("Frequency", 2, line=3.5, cex=1.5)
par(new=T)
hist(mod2, main="",xlim=c(0,100), ylim=c(0,60), xlab="", breaks=35, col="gray77", yaxt="n", xaxt="n", ylab="")
legend("topright", legend=c("Model 1","Model 2"), fill=c("white","gray77"), cex=2.5)

par(mai=c(1,.1,.4,.1)) # BLTR
dotchart(c(1.184227,0.278557), labels=c("Model 2","Model 1"), xlim=c(-.05,1.5), pch=19, cex=1.5)
segments(coef(res1)[1]+0.00833*1.96, 2, coef(res1)[1]-0.00833*1.96, 2, lwd=3)
segments(coef(res1)[2]+0.01373*1.96, 1, coef(res1)[2]-0.01373*1.96, 1, lwd=3)

#dev.off()


##############
## FIGURE 3 ##

## Apply Slapin's positioning and error to a 10 actor model ###


# NBS Function for 10 actors #
# s is sq, positions of actors m1 to m10 represented as a1 to a10. Comission position is a11. We seek to maximize c1* ... * c10 * com (combined utility) by varying position x

nb <- function(s,x,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11){
  c1 <- ((-1*(a1-x))^2)-((-1*(a1-s))^2)
  c2 <- ((-1*(a2-x))^2)-((-1*(a2-s))^2)
  c3 <- ((-1*(a3-x))^2)-((-1*(a3-s))^2)
  c4 <- ((-1*(a4-x))^2)-((-1*(a4-s))^2)
  c5 <- ((-1*(a5-x))^2)-((-1*(a5-s))^2)
  c6 <- ((-1*(a6-x))^2)-((-1*(a6-s))^2)
  c7 <- ((-1*(a7-x))^2)-((-1*(a7-s))^2)
  c8 <- ((-1*(a8-x))^2)-((-1*(a8-s))^2)
  c9 <- ((-1*(a9-x))^2)-((-1*(a9-s))^2)
  c10 <- ((-1*(a10-x))^2)-((-1*(a10-s))^2)
  com <- ((-1*(a11-x))^2)-((-1*(a11-s))^2)
  c1*c2*c3*c4*c5*c6*c7*c8*c9*c10*com
}

# Basics
nissues <- 150
sd <- c(2,12)

# Create results objects
mean.abs.diff.RR.JS.10actors <- matrix(NA,2,1000)
mean.abs.diff.NBS.JS.10actors <- matrix(NA,2,1000)
mean.abs.diff.CM.JS.10actors <- matrix(NA,2,1000)

# Correlations in model predictions
cor.RR.NBS.10actors <- matrix(NA,2,1000)
cor.RR.CM.10actors <- matrix(NA,2,1000)


for (k in 1:length(sd)){

    for (j in 1:1000){
        # Create true positions and outcome
        m1 <- round(runif(nissues, min=0, max=60))
        m2 <- round(runif(nissues, min=0, max=60))
        m3 <- round(runif(nissues, min=0, max=60))
        m4 <- round(runif(nissues, min=0, max=60))
        m5 <- round(runif(nissues, min=0, max=60))
        m6 <- round(runif(nissues, min=0, max=60))
        m7 <- round(runif(nissues, min=0, max=60))
        m8 <- round(runif(nissues, min=0, max=60))
        m9 <- round(runif(nissues, min=0, max=60))
        m10 <- round(runif(nissues, min=0, max=60))
        x2 <- round(runif(nissues, min=40, max=100))
    
        issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
        true.pivot <- rep(NA,150)
        for (i in 1:150){
            true.pivot[i] <- sort(issue.positions[i,1:10])[4]
        }
        x1 <- true.pivot
        x2 <- ifelse(x1==x2,x2+1,x2) # Ensure the two positions do not agree
        outtrue <- ifelse((2*x1)<x2,2*x1,x2)

        # Create observed positions, outcomes, and sq (our measurements)
        m1.obs <- round(m1+rnorm(nissues,0,sd[k]),0)
        m2.obs <- round(m2+rnorm(nissues,0,sd[k]),0)
        m3.obs <- round(m3+rnorm(nissues,0,sd[k]),0)
        m4.obs <- round(m4+rnorm(nissues,0,sd[k]),0)
        m5.obs <- round(m5+rnorm(nissues,0,sd[k]),0)
        m6.obs <- round(m6+rnorm(nissues,0,sd[k]),0)
        m7.obs <- round(m7+rnorm(nissues,0,sd[k]),0)
        m8.obs <- round(m8+rnorm(nissues,0,sd[k]),0)
        m9.obs <- round(m9+rnorm(nissues,0,sd[k]),0)
        m10.obs <- round(m10+rnorm(nissues,0,sd[k]),0)
        x2.obs <- round(x2+rnorm(nissues,0,sd[k]),0)
        out.obs <- round(outtrue+rnorm(nissues,0,sd[k]),0)
        sq.obs  <- round(rnorm(nissues,0,sd[k]),0)

        # Rescale to place on 0-100 scale
        absmax <- max(apply(cbind(m1.obs,m2.obs,m3.obs,m4.obs,m5.obs,m6.obs,m7.obs,m8.obs,m9.obs,m10.obs,x2.obs,out.obs,sq.obs),1,max))
        absmin<- min(apply(cbind(m1.obs,m2.obs,m3.obs,m4.obs,m5.obs,m6.obs,m7.obs,m8.obs,m9.obs,m10.obs,x2.obs,out.obs,sq.obs),1,min))
        range <- absmax-absmin
        m1.obs <- round(((m1.obs+abs(absmin))/range)*100,0)
        m2.obs <- round(((m2.obs+abs(absmin))/range)*100,0)
        m3.obs <- round(((m3.obs+abs(absmin))/range)*100,0)
        m4.obs <- round(((m4.obs+abs(absmin))/range)*100,0)
        m5.obs <- round(((m5.obs+abs(absmin))/range)*100,0)
        m6.obs <- round(((m6.obs+abs(absmin))/range)*100,0)
        m7.obs <- round(((m7.obs+abs(absmin))/range)*100,0)
        m8.obs <- round(((m8.obs+abs(absmin))/range)*100,0)
        m9.obs <- round(((m9.obs+abs(absmin))/range)*100,0)
        m10.obs <- round(((m10.obs+abs(absmin))/range)*100,0)
        x2.obs <- round(((x2.obs+abs(absmin))/range)*100,0)

        out.obs <- round(((out.obs+abs(absmin))/range)*100,0)
        sq.obs <- round(((sq.obs+abs(absmin))/range)*100,0)

    
        # R-R model applied to observed data
        issue.positions.obs <- cbind(m1.obs,m2.obs,m3.obs,m4.obs,m5.obs,m6.obs,m7.obs,m8.obs,m9.obs,m10.obs)
        pivot.obs <- rep(NA,150)
        for (i in 1:150){
            pivot.obs[i] <- sort(issue.positions.obs[i,1:10])[4]
        }
        x1.obs <- pivot.obs
        x2.obs <- ifelse(x1.obs==x2.obs,x2.obs+1,x2.obs) # Need to make sure that the two positions do not agree

        # For simplicity, ensure sq is never to left of both observed positions.
        sq.obs[sq.obs>x1.obs & sq.obs>x2.obs] <- x2.obs[sq.obs>x1.obs & sq.obs>x2.obs] - 1


        outrrobs <- rep(NA,nissues) # Holding bin for R-R observations
        # create scenarios --- SQ inbetween observed positions
        outrrobs[x1.obs<=sq.obs & sq.obs<=x2.obs | x2.obs<=sq.obs & sq.obs<=x1.obs]<- sq.obs[x1.obs<=sq.obs & sq.obs<=x2.obs | x2.obs<=sq.obs & sq.obs<=x1.obs]
        # create scenarios --- SQ to the right of both positions
        outrrobs[x2.obs>sq.obs & x1.obs>sq.obs & 2*(x1.obs-sq.obs)>=(x2.obs-sq.obs)] <- x2.obs[x2.obs>sq.obs & x1.obs>sq.obs & 2*(x1.obs-sq.obs)>=(x2.obs-sq.obs)]
        outrrobs[x2.obs>sq.obs & x1.obs>sq.obs & 2*(x1.obs-sq.obs)<(x2.obs-sq.obs)]<-2*x1.obs[x2.obs>sq.obs & x1.obs>sq.obs & 2*(x1.obs-sq.obs)<(x2.obs-sq.obs)]
        # create scenarios --- SQ to the left of both positions

        # Calculate mean(abs(difference)) - R-R
        mean.abs.diff.RR.JS.10actors[k,j] <- mean(abs(outtrue-outrrobs))

        # NB applied to simulated data
        outnbobs <- rep(NA,nissues)
        for (i in 1:nissues){
            res <- optimize(nb, lower=0, upper=100, a1=m1.obs[i], a2=m2.obs[i], a3=m3.obs[i], a4=m4.obs[i], a5=m5.obs[i], a6=m6.obs[i], a7=m7.obs[i], a8=m8.obs[i], a9=m9.obs[i], a10=m10.obs[i], a11=x2.obs[i], s=sq.obs[i], maximum=TRUE)
            outnbobs[i] <- res$maximum	
        }
    
        # Calculate mean(abs(difference)) - NBS
        mean.abs.diff.NBS.JS.10actors[k,j] <- mean(abs(outtrue-outnbobs))

        # CM applied to simulated data
        outcmobs <- rep(NA, nissues)
        outcmobs <- (m1.obs+m2.obs+m3.obs+m4.obs+m5.obs+m6.obs+m7.obs+m8.obs+m9.obs+m10.obs+x2.obs)/11

        # Calculate mean(abs(difference)) - CM
        mean.abs.diff.CM.JS.10actors[k,j] <- mean(abs(outtrue-outcmobs))
    
        # Check correlations between predictions
        cor.RR.NBS.10actors[k,j] <- cor(outrrobs,outnbobs)
        cor.RR.CM.10actors[k,j] <- cor(outrrobs,outcmobs)
    }
}




#jpeg("Figures/Figure-3.jpg", width=900, height=600)
par(mfrow=c(1,2))
boxplot(mean.abs.diff.RR.JS.10actors[1,], mean.abs.diff.NBS.JS.10actors[1,], mean.abs.diff.CM.JS.10actors[1,], ylim=c(0,40), names=c("RR","NBS","CM"), main="", cex.axis=2.2, xaxt="n")
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Measurement Error: sd = 2", 3, cex=2, line=1)
boxplot(mean.abs.diff.RR.JS.10actors[2,], mean.abs.diff.NBS.JS.10actors[2,], mean.abs.diff.CM.JS.10actors[2,], ylim=c(0,40), names=c("RR","NBS","CM"), main="", cex.axis=2.2, xaxt="n")
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Measurement Error: sd = 12", 3, cex=2, line=1)
#dev.off()





##############
## FIGURE 4 ##


# NBS Function #
# s is sq, positions of actors m1 to m10 represented as a1 to a10. Comission position is a11. We seek to maximize c1* ... * c10 * com (combined utility) by varying position x

nb <- function(s,x,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11){
  c1 <- ((-1*(a1-x))^2)-((-1*(a1-s))^2)
  c2 <- ((-1*(a2-x))^2)-((-1*(a2-s))^2)
  c3 <- ((-1*(a3-x))^2)-((-1*(a3-s))^2)
  c4 <- ((-1*(a4-x))^2)-((-1*(a4-s))^2)
  c5 <- ((-1*(a5-x))^2)-((-1*(a5-s))^2)
  c6 <- ((-1*(a6-x))^2)-((-1*(a6-s))^2)
  c7 <- ((-1*(a7-x))^2)-((-1*(a7-s))^2)
  c8 <- ((-1*(a8-x))^2)-((-1*(a8-s))^2)
  c9 <- ((-1*(a9-x))^2)-((-1*(a9-s))^2)
  c10 <- ((-1*(a10-x))^2)-((-1*(a10-s))^2)
  com <- ((-1*(a11-x))^2)-((-1*(a11-s))^2)
  c1*c2*c3*c4*c5*c6*c7*c8*c9*c10*com
}

############################
# Actors Misallocated: 10% #

# Objects to capture mean absolute differences
mean.abs.diff.RR.1actor <- rep(NA,1000)
mean.abs.diff.NB.1actor <- rep(NA,1000)
mean.abs.diff.CM.1actor <- rep(NA,1000)


for (j in 1:1000){

# Calculate outtrue (true R-R outcomes, no uncertainty)
  nissues <-150 # number of issues
  m1 <- sample(c(0,50,100), 150, replace=T)
  m2 <- sample(c(0,50,100), 150, replace=T)
  m3 <- sample(c(0,50,100), 150, replace=T)
  m4 <- sample(c(0,50,100), 150, replace=T)
  m5 <- sample(c(0,50,100), 150, replace=T)
  m6 <- sample(c(0,50,100), 150, replace=T)
  m7 <- sample(c(0,50,100), 150, replace=T)
  m8 <- sample(c(0,50,100), 150, replace=T)
  m9 <- sample(c(0,50,100), 150, replace=T)
  m10 <- sample(c(0,50,100), 150, replace=T)
  x2 <-  sample(c(50,100), 150, replace=T)
  issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  true.pivot <- rep(NA,150)
  for (i in 1:150){
    true.pivot[i] <- sort(issue.positions[i,1:10])[4]
  }
  x1 <- true.pivot
  outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Add uncertainty to the positioning of actors: 1 out of the 10 allocated to wrong group (Remake m1)
  m1.redo <- rep(NA,150)
  
  for (i in 1:150){
    if(m1[i]==0) {
      m1.redo[i] <- sample(c(50,100),1)
    } else {
      if(m1[i]==50){
        m1.redo[i] <- sample(c(0,100),1)
      } else {
        if(m1[i]==100){
          m1.redo[i] <- sample(c(0,50),1)
        }
      }
    }
  }

# R-R applied to misallocated data (X2 is agenda setter)
# Re-calculate the pivot given these ten actors' positions
  issue.pos.redo <- cbind(m1.redo,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  pivot.redo <- rep(NA,150)
  for (i in 1:150){
    pivot.redo[i] <- sort(issue.pos.redo[i,1:10])[4]
  }
  x1.redo <- pivot.redo

  outtrue.redo <- ifelse((2*x1.redo)<x2,2*x1.redo,x2)

# NBS applied to misallocated data
  outnb.redo <- rep(NA,nissues)
  for (i in 1:nissues){
    res <- optimize(nb, lower=0, upper=100, a1=m1.redo[i], a2=m2[i], a3=m3[i], a4=m4[i], a5=m5[i], a6=m6[i], a7=m7[i], a8=m8[i], a9=m9[i], a10=m10[i], a11=x2[i], s=0, maximum=TRUE)
    outnb.redo[i] <- res$maximum	
  }

# CM applied to misallocated data
  outcm.redo <- (m1.redo+m2+m3+m4+m5+m6+m7+m8+m9+m10+x2)/11

# Calculate mean(abs(difference))
  mean.abs.diff.RR.1actor[j] <- mean(abs(outtrue-outtrue.redo))
  mean.abs.diff.NB.1actor[j] <- mean(abs(outtrue-outnb.redo))
  mean.abs.diff.CM.1actor[j] <- mean(abs(outtrue-outcm.redo))
}

boxplot(mean.abs.diff.RR.1actor, mean.abs.diff.NB.1actor,mean.abs.diff.CM.1actor, ylim=c(0,60), names=c("R-R","NBS","CM"), main="Actors Misallocated: 10%")


############################
# Actors Misallocated: 20% #

# Objects to capture mean absolute differences
mean.abs.diff.RR.2actor <- rep(NA,1000)
mean.abs.diff.NB.2actor <- rep(NA,1000)
mean.abs.diff.CM.2actor <- rep(NA,1000)

for (j in 1:1000){

# Calculate outtrue (true R-R outcomes, no uncertainty)
  nissues <-150
  m1 <- sample(c(0,50,100), 150, replace=T)
  m2 <- sample(c(0,50,100), 150, replace=T)
  m3 <- sample(c(0,50,100), 150, replace=T)
  m4 <- sample(c(0,50,100), 150, replace=T)
  m5 <- sample(c(0,50,100), 150, replace=T)
  m6 <- sample(c(0,50,100), 150, replace=T)
  m7 <- sample(c(0,50,100), 150, replace=T)
  m8 <- sample(c(0,50,100), 150, replace=T)
  m9 <- sample(c(0,50,100), 150, replace=T)
  m10 <- sample(c(0,50,100), 150, replace=T)
  x2 <-  sample(c(50,100), 150, replace=T)
  issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  true.pivot <- rep(NA,150)
  for (i in 1:150){
    true.pivot[i] <- sort(issue.positions[i,1:10])[4]
  }
  x1 <- true.pivot
  outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Add uncertainty to the positioning of actors: 2 out of the 10 allocated to wrong group (Remake m1 and m2)

  m1.redo <- rep(NA,150)
  for (i in 1:150){
    if(m1[i]==0) {
      m1.redo[i] <- sample(c(50,100),1)
    } else {
      if(m1[i]==50){
        m1.redo[i] <- sample(c(0,100),1)
      } else {
        if(m1[i]==100){
          m1.redo[i] <- sample(c(0,50),1)
        }
      }
    }
  }

  m2.redo <- rep(NA,150)
  for (i in 1:150){
    if(m2[i]==0) {
      m2.redo[i] <- sample(c(50,100),1)
    } else {
      if(m2[i]==50){
        m2.redo[i] <- sample(c(0,100),1)
      } else {
        if(m2[i]==100){
          m2.redo[i] <- sample(c(0,50),1)
        }
      }
    }
  }

# R-R applied to misallocated data (X2 is agenda setter)
# Re-calculate the pivot given these ten actors' positions
  issue.pos.redo.2 <- cbind(m1.redo,m2.redo,m3,m4,m5,m6,m7,m8,m9,m10)
  pivot.redo.2 <- rep(NA,150)
  for (i in 1:150){
    pivot.redo.2[i] <- sort(issue.pos.redo.2[i,1:10])[4]
  }
  x1.redo.2 <- pivot.redo.2
  outtrue.redo.2 <- ifelse((2*x1.redo.2)<x2,2*x1.redo.2,x2)

# NBS applied to misallocated data
  outnb.redo.2 <- rep(NA,nissues)
  for (i in 1:nissues){
    res <- optimize(nb, lower=0, upper=100, a1=m1.redo[i], a2=m2.redo[i], a3=m3[i], a4=m4[i], a5=m5[i], a6=m6[i], a7=m7[i], a8=m8[i], a9=m9[i], a10=m10[i], a11=x2[i], s=0, maximum=TRUE)
    outnb.redo.2[i] <- res$maximum	
  }

# CM applied to misallocated data
  outcm.redo.2 <- (m1.redo+m2.redo+m3+m4+m5+m6+m7+m8+m9+m10+x2)/11

# Calculate mean(abs(difference))
  mean.abs.diff.RR.2actor[j] <- mean(abs(outtrue-outtrue.redo.2))
  mean.abs.diff.NB.2actor[j] <- mean(abs(outtrue-outnb.redo.2))
  mean.abs.diff.CM.2actor[j] <- mean(abs(outtrue-outcm.redo.2))
}


boxplot(mean.abs.diff.RR.2actor, mean.abs.diff.NB.2actor,mean.abs.diff.CM.2actor, ylim=c(0,60), names=c("R-R","NBS","CM"), main="Actors Misallocated: 20%")



########################################
# Error in Relative Positioning, sd=10 #
# Clusters correctly identified, but incorrectly placed on dimension. Score of the middle group varies by sd = 10.

mean.abs.diff.RR.sd10 <- rep(NA,1000)
mean.abs.diff.NB.sd10 <- rep(NA,1000)
mean.abs.diff.CM.sd10 <- rep(NA,1000)

for (j in 1:1000){

# Calculate outtrue
  nissues <-150
  m1 <- sample(c(0,50,100), 150, replace=T)
  m2 <- sample(c(0,50,100), 150, replace=T)
  m3 <- sample(c(0,50,100), 150, replace=T)
  m4 <- sample(c(0,50,100), 150, replace=T)
  m5 <- sample(c(0,50,100), 150, replace=T)
  m6 <- sample(c(0,50,100), 150, replace=T)
  m7 <- sample(c(0,50,100), 150, replace=T)
  m8 <- sample(c(0,50,100), 150, replace=T)
  m9 <- sample(c(0,50,100), 150, replace=T)
  m10 <- sample(c(0,50,100), 150, replace=T)
  x2 <-  sample(c(50,100), 150, replace=T)

# Create True Romer-Rosenthal Outcome (x2 is agenda setter)
  issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  true.pivot <- rep(NA,150)
  for (i in 1:150){
    true.pivot[i] <- sort(issue.positions[i,1:10])[4]
  }
  x1 <- true.pivot
  outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Add error to all positions in middle cluster
  SD.10 <- 50+rnorm(150,0,10)
  issue.pos <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,x2,outtrue)
  issue.pos.redo.10 <- matrix(NA,150,12)
  for (i in 1:150){
    issue.pos.redo.10[i,] <- ifelse(issue.pos[i,]==50, SD.10[i],issue.pos[i,])
  }

# Scale any middle clusters that climbed above 100 or below 0 to the fixed endpoints.
# Replace values above 100 with 100
  test.a <- matrix(NA,150,12)
  for (i in 1:150){
    test.a[i,] <- ifelse(issue.pos.redo.10[i,]>100, 100, issue.pos.redo.10[i,])
  }

# Replace values below 0 with 0
  test.b <- matrix(NA,150,12)
  for (i in 1:150){
    test.b[i,] <- ifelse(test.a[i,]<0, 0, test.a[i,])
  }
  issue.pos.redo.10 <- test.b

# R-R applied to misallocated data (x2 is agenda setter)
# Re-calculate the pivot given these ten actors' positions
  pivot.redo.10 <- rep(NA,150)
  for (i in 1:150){
    pivot.redo.10[i] <- sort(issue.pos.redo.10[i,1:10])[4]
  }
  x1.redo.10 <- pivot.redo.10
  x2.redo.10 <- issue.pos.redo.10[,11]
  outtrue.redo.10 <- ifelse((2*x1.redo.10)<x2.redo.10,2*x1.redo.10,x2.redo.10)
  
# NBS applied to misidentified cluster positions
  outnb.redo.10 <- rep(NA,nissues)
  for (i in 1:nissues){
    res <- optimize(nb, lower=0, upper=100, a1=issue.pos.redo.10[i,1], a2=issue.pos.redo.10[i,2], a3=issue.pos.redo.10[i,3], a4=issue.pos.redo.10[i,4], a5=issue.pos.redo.10[i,5], a6=issue.pos.redo.10[i,6], a7=issue.pos.redo.10[i,7], a8=issue.pos.redo.10[i,8], a9=issue.pos.redo.10[i,9], a10=issue.pos.redo.10[i,10], a11=issue.pos.redo.10[i,11], s=0, maximum=TRUE)
    outnb.redo.10[i] <- res$maximum	
  }

# CM applied to misidentified cluster positions
  outcm.redo.10 <- rowSums(issue.pos.redo.10[,1:11])/11

# outtrue modified by same error because no uncertainty between outcome position and middle cluster position
  outtrue.orig.10 <- issue.pos.redo.10[,12]


# Calculate mean(abs(difference))
  mean.abs.diff.RR.sd10[j] <- mean(abs(outtrue.orig.10-outtrue.redo.10))
  mean.abs.diff.NB.sd10[j] <- mean(abs(outtrue.orig.10-outnb.redo.10))
  mean.abs.diff.CM.sd10[j] <- mean(abs(outtrue.orig.10-outcm.redo.10))
}


boxplot(mean.abs.diff.RR.sd10, mean.abs.diff.NB.sd10,mean.abs.diff.CM.sd10, ylim=c(0,60), names=c("R-R","NB","CM"), main="Error in Relative Positioning, sd=10")




########################################
# Error in Relative Positioning, sd=20 #
# Clusters correctly identified, but incorrectly placed on dimension. Score of the middle group varies by sd = 20.

mean.abs.diff.RR.sd20 <- rep(NA,1000)
mean.abs.diff.NB.sd20 <- rep(NA,1000)
mean.abs.diff.CM.sd20 <- rep(NA,1000)

for (j in 1:1000){

# Calculate outtrue
  nissues <-150
  m1 <- sample(c(0,50,100), 150, replace=T)
  m2 <- sample(c(0,50,100), 150, replace=T)
  m3 <- sample(c(0,50,100), 150, replace=T)
  m4 <- sample(c(0,50,100), 150, replace=T)
  m5 <- sample(c(0,50,100), 150, replace=T)
  m6 <- sample(c(0,50,100), 150, replace=T)
  m7 <- sample(c(0,50,100), 150, replace=T)
  m8 <- sample(c(0,50,100), 150, replace=T)
  m9 <- sample(c(0,50,100), 150, replace=T)
  m10 <- sample(c(0,50,100), 150, replace=T)
  x2 <-  sample(c(50,100), 150, replace=T)

# Create True Romer-Rosenthal Outcome (x2 is agenda setter)
  issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  true.pivot <- rep(NA,150)
  for (i in 1:150){
    true.pivot[i] <- sort(issue.positions[i,1:10])[4]
  }
  x1 <- true.pivot
  outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Add error to all positions in middle cluster
  SD.20 <- 50+rnorm(150,0,20)
  issue.pos <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,x2,outtrue)
  issue.pos.redo.20 <- matrix(NA,150,12)
  for (i in 1:150){
    issue.pos.redo.20[i,] <- ifelse(issue.pos[i,]==50, SD.20[i],issue.pos[i,])
  }
  
# Scale any middle clusters that climbed above 100 or below 0 to the fixed endpoints.
# Replace values above 100 with 100
  test.a <- matrix(NA,150,12)
  for (i in 1:150){
    test.a[i,] <- ifelse(issue.pos.redo.20[i,]>100, 100, issue.pos.redo.20[i,])
  }

# Replace values below 0 with 0
  test.b <- matrix(NA,150,12)
  for (i in 1:150){
    test.b[i,] <- ifelse(test.a[i,]<0, 0, test.a[i,])
  }
  issue.pos.redo.20 <- test.b

# R-R applied to misallocated data (x2 is agenda setter)
# Re-calculate the pivot given these ten actors' positions
  pivot.redo.20 <- rep(NA,150)

  for (i in 1:150){
    pivot.redo.20[i] <- sort(issue.pos.redo.20[i,1:10])[4]
  }

  x1.redo.20 <- pivot.redo.20
  x2.redo.20 <- issue.pos.redo.20[,11]
  outtrue.redo.20 <- ifelse((2*x1.redo.20)<x2.redo.20,2*x1.redo.20,x2.redo.20)

# NBS applied to misidentified cluster positions
  outnb.redo.20 <- rep(NA,nissues)
  for (i in 1:nissues){
    res <- optimize(nb, lower=0, upper=100, a1=issue.pos.redo.20[i,1], a2=issue.pos.redo.20[i,2], a3=issue.pos.redo.20[i,3], a4=issue.pos.redo.20[i,4], a5=issue.pos.redo.20[i,5], a6=issue.pos.redo.20[i,6], a7=issue.pos.redo.20[i,7], a8=issue.pos.redo.20[i,8], a9=issue.pos.redo.20[i,9], a10=issue.pos.redo.20[i,10], a11=issue.pos.redo.20[i,11], s=0, maximum=TRUE)
    outnb.redo.20[i] <- res$maximum	
  }
  
# CM applied to misidentified cluster positions
  outcm.redo.20 <- rowSums(issue.pos.redo.20[,1:11])/11
  
# Outtrue modified by same error because no uncertainty between outcome position and middle cluster position
  outtrue.orig.20 <- issue.pos.redo.20[,12]
  
# Calculate mean(abs(difference))
  mean.abs.diff.RR.sd20[j] <- mean(abs(outtrue.orig.20-outtrue.redo.20))
  mean.abs.diff.NB.sd20[j] <- mean(abs(outtrue.orig.20-outnb.redo.20))
  mean.abs.diff.CM.sd20[j] <- mean(abs(outtrue.orig.20-outcm.redo.20))
}

boxplot(mean.abs.diff.RR.sd20, mean.abs.diff.NB.sd20,mean.abs.diff.CM.sd20, ylim=c(0,60), names=c("R-R","NBS","CM"), main="Error in Relative Positioning, sd=20")



###################################################################
# Actors Misallocated: 20% & Error in Relative Positioning, sd=20 #

mean.abs.diff.RR.combo <- rep(NA,1000)
mean.abs.diff.NB.combo <- rep(NA,1000)
mean.abs.diff.CM.combo <- rep(NA,1000)

for (j in 1:1000){

# Calculate outtrue
  nissues <-150
  m1 <- sample(c(0,50,100), 150, replace=T)
  m2 <- sample(c(0,50,100), 150, replace=T)
  m3 <- sample(c(0,50,100), 150, replace=T)
  m4 <- sample(c(0,50,100), 150, replace=T)
  m5 <- sample(c(0,50,100), 150, replace=T)
  m6 <- sample(c(0,50,100), 150, replace=T)
  m7 <- sample(c(0,50,100), 150, replace=T)
  m8 <- sample(c(0,50,100), 150, replace=T)
  m9 <- sample(c(0,50,100), 150, replace=T)
  m10 <- sample(c(0,50,100), 150, replace=T)
  x2 <-  sample(c(50,100), 150, replace=T)

# Create True Romer-Rosenthal Outcome, X2 is agenda setter
  issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  true.pivot <- rep(NA,150)
  for (i in 1:150){
    true.pivot[i] <- sort(issue.positions[i,1:10])[4]
  }
  x1 <- true.pivot
  outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Positioning of actors: 2 out of the 10 allocated to wrong group (Remake m1 and m2)
  m1.redo <- rep(NA,150)
  for (i in 1:150){
    if(m1[i]==0) {
      m1.redo[i] <- sample(c(50,100),1)
    } else {
      if(m1[i]==50){
        m1.redo[i] <- sample(c(0,100),1)
      } else {
        if(m1[i]==100){
          m1.redo[i] <- sample(c(0,50),1)
        }
      }
    }
  }

  m2.redo <- rep(NA,150)
  for (i in 1:150){
    if(m2[i]==0) {
      m2.redo[i] <- sample(c(50,100),1)
    } else {
      if(m2[i]==50){
        m2.redo[i] <- sample(c(0,100),1)
      } else {
        if(m2[i]==100){
          m2.redo[i] <- sample(c(0,50),1)
        }
      }
    }
  }

# Add error to all positions in middle cluster
  SD.20 <- 50+rnorm(150,0,20)
  issue.pos.combo <- cbind(m1.redo,m2.redo,m3,m4,m5,m6,m7,m8,m9,m10,x2,outtrue)
  issue.pos.combo.20 <- matrix(NA,150,12)
  for (i in 1:150){
    issue.pos.combo.20[i,] <- ifelse(issue.pos.combo[i,]==50, SD.20[i], issue.pos.combo[i,])
  }

# Scale any middle clusters that climbed above 100 or below 0 to the fixed endpoints.
# Replace values above 100 with 100
  test.a <- matrix(NA,150,12)
  for (i in 1:150){
    test.a[i,] <- ifelse(issue.pos.combo.20[i,]>100, 100, issue.pos.combo.20[i,])
  }

# Replace values below 0 with 0
  test.b <- matrix(NA,150,12)
  for (i in 1:150){
    test.b[i,] <- ifelse(test.a[i,]<0, 0, test.a[i,])
  }
  issue.pos.combo.20 <- test.b

# R-R applied to misallocated data (x2 is agenda setter)
# Re-calculate the pivot given these ten actors' positions
  pivot.combo.20 <- rep(NA,150)
  for (i in 1:150){
    pivot.combo.20[i] <- sort(issue.pos.combo.20[i,1:10])[4]
  }
  x1.combo.20 <- pivot.combo.20
  x2.combo.20 <- issue.pos.combo.20[,11]

# R-R applied to misidentified cluster positions and 2 misallocated actors
  outtrue.combo.20 <- ifelse((2*x1.combo.20)<x2.combo.20, 2*x1.combo.20, x2.combo.20)

# NBS applied to misidentified cluster positions and 2 misallocated actors
  outnb.combo.20 <- rep(NA,nissues)
  for (i in 1:nissues){
    res <- optimize(nb, lower=0, upper=100, a1=issue.pos.combo.20[i,1], a2=issue.pos.combo.20[i,2], a3=issue.pos.combo.20[i,3], a4=issue.pos.combo.20[i,4], a5=issue.pos.combo.20[i,5], a6=issue.pos.combo.20[i,6], a7=issue.pos.combo.20[i,7], a8=issue.pos.combo.20[i,8], a9=issue.pos.combo.20[i,9], a10=issue.pos.combo.20[i,10], a11=issue.pos.combo.20[i,11], s=0, maximum=TRUE)
    outnb.combo.20[i] <- res$maximum	
  }

# CM applied to misidentified cluster positions and 2 misallocated actors
  outcm.combo.20 <- rowSums(issue.pos.combo.20[,1:11])/11

# Outtrue modified
  outtrue.orig.20 <- issue.pos.combo.20[,12]

# Calculate mean(abs(difference))
  mean.abs.diff.RR.combo[j] <- mean(abs(outtrue.orig.20-outtrue.combo.20))
  mean.abs.diff.NB.combo[j] <- mean(abs(outtrue.orig.20-outnb.combo.20))
  mean.abs.diff.CM.combo[j] <- mean(abs(outtrue.orig.20-outcm.combo.20))
}

boxplot(mean.abs.diff.RR.combo, mean.abs.diff.NB.combo, mean.abs.diff.CM.combo, ylim=c(0,60), names=c("R-R","NBs","CM"), main="Error in Allocation (20%) & Positioning (sd = 20)")



###################################################################
# Random error applied to any 1 per row of: states, commission, SQ and outcome #

mean.abs.diff.RR.random <- rep(NA,1000)
mean.abs.diff.NB.random <- rep(NA,1000)
mean.abs.diff.CM.random <- rep(NA,1000)

for (j in 1:1000){

# Calculate outtrue (true RR outcomes, no uncertainty)
  nissues <-150
  m1 <- sample(c(0,50,100), 150, replace=T)
  m2 <- sample(c(0,50,100), 150, replace=T)
  m3 <- sample(c(0,50,100), 150, replace=T)
  m4 <- sample(c(0,50,100), 150, replace=T)
  m5 <- sample(c(0,50,100), 150, replace=T)
  m6 <- sample(c(0,50,100), 150, replace=T)
  m7 <- sample(c(0,50,100), 150, replace=T)
  m8 <- sample(c(0,50,100), 150, replace=T)
  m9 <- sample(c(0,50,100), 150, replace=T)
  m10 <- sample(c(0,50,100), 150, replace=T)
  x2 <-  sample(c(50,100), 150, replace=T)
  issue.positions <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)
  true.pivot <- rep(NA,150)
  for (i in 1:150){
    true.pivot[i] <- sort(issue.positions[i,1:10])[4]
  }
  x1 <- true.pivot
  outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Introduce 1 piece of random error in allocation per issue (row).
# Randomly select 1 cell per row and replace with unchosen option (0, 50, 100)
  orig.mat <- cbind(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,x2,outtrue)
  error.mat <- orig.mat
  sample1 <- sample(1:12,150, replace=T)

  for (i in 1:150){
    if(orig.mat[i,sample1[i]]==0) {
      error.mat[i,sample1[i]] <- sample(c(50,100),1)
    } else {
      if(orig.mat[i,sample1[i]]==50) {
        error.mat[i,sample1[i]] <- sample(c(0,100),1)
      } else {
        if(orig.mat[i,sample1[i]]==100) {
          error.mat[i,sample1[i]] <- sample(c(0,50),1)
        }
      }
    }
  }

# Re-calculate the pivot given these ten actors' positions
  pivot.redo <- rep(NA,150)
  for (i in 1:150){
    pivot.redo[i] <- sort(error.mat[i,1:10])[4]
  }
  x1.redo <- pivot.redo
  x2.redo <- error.mat[,11]
  
# RR applied to misallocated data
  outtrue.redo <- ifelse((2*x1.redo)<x2.redo,2*x1.redo,x2.redo)

# NB applied to misallocated data
  outnb.redo <- rep(NA,nissues)
  for (i in 1:nissues){
    res <- optimize(nb, lower=0, upper=100, a1=error.mat[i,1], a2=error.mat[i,2], a3=error.mat[i,3], a4=error.mat[i,4], a5=error.mat[i,5], a6=error.mat[i,6], a7=error.mat[i,7], a8=error.mat[i,8], a9=error.mat[i,9], a10=error.mat[i,10], a11=error.mat[i,11], s=0, maximum=TRUE)
    outnb.redo[i] <- res$maximum	
  }

# CM applied to data with random error
  outcm.redo <-rowSums(error.mat[,1:11])/11
  
# Outtrue modified
  outtrue.error <- error.mat[,12]
  
# Calculate mean(abs(difference))
  mean.abs.diff.RR.random[j] <- mean(abs(outtrue.error-outtrue.redo))
  mean.abs.diff.NB.random[j] <- mean(abs(outtrue.error-outnb.redo))
  mean.abs.diff.CM.random[j] <- mean(abs(outtrue.error-outcm.redo))
}


boxplot(mean.abs.diff.RR.random, mean.abs.diff.NB.random, mean.abs.diff.CM.random, ylim=c(0,60), names=c("R-R","NBS","CM"), main="Misallocated Actor, SQ or Outcome: Random Selection", cex=1.5, cex.axis=1.5, cex.main=1.5)



#jpeg("Figures/Figure-4.jpg", width=768, height=1024, quality=100)
nf <- layout(matrix(c(1,1,2,2,3,3,4,4,5,5,6,6), 3,4, byrow=T))
#layout.show(nf)

boxplot(mean.abs.diff.RR.1actor, mean.abs.diff.NB.1actor,mean.abs.diff.CM.1actor, ylim=c(0,60), names=c("R-R","NBS","CM"), main="", cex.axis=2.2, xaxt="n")
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Actors Misallocated: 10%", 3, cex=1.5, line=1)

boxplot(mean.abs.diff.RR.2actor, mean.abs.diff.NB.2actor,mean.abs.diff.CM.2actor, ylim=c(0,60), names=c("R-R","NBS","CM"), main="", cex.axis=2.2, xaxt="n")
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Actors Misallocated: 20%", 3, cex=1.5, line=1)

boxplot(mean.abs.diff.RR.sd10, mean.abs.diff.NB.sd10,mean.abs.diff.CM.sd10, ylim=c(0,60), names=c("R-R","NBS","CM"), main="", cex.axis=2.2, xaxt="n")
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Error in Relative Positioning: sd = 10", 3, cex=1.5, line=1)

boxplot(mean.abs.diff.RR.sd20, mean.abs.diff.NB.sd20,mean.abs.diff.CM.sd20, ylim=c(0,60), names=c("R-R","NBS","CM"), main="", cex.axis=2.2, xaxt="n")
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Error in Relative Positioning: sd = 20", 3, cex=1.5, line=1)

boxplot(mean.abs.diff.RR.combo, mean.abs.diff.NB.combo,mean.abs.diff.CM.combo, ylim=c(0,60), names=c("R-R","NBS","CM"), main="", cex.axis=2.2, xaxt="n", mai=c(1,2,1,2))
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Error in Allocation (20%) & Positioning (sd = 20)", 3, cex=1.5, line=1)

boxplot(mean.abs.diff.RR.random, mean.abs.diff.NB.random, mean.abs.diff.CM.random, ylim=c(0,60), names=c("R-R","NBS","CM"), main="", cex.axis=2.2, xaxt="n", mai=c(1,2,1,2))
axis(1, at=1:3, tick=F, labels=c("R-R","NBS","CM"), cex.axis=2.2, line=1)
mtext("Misallocated Actor, SQ or Outcome", 3, cex=1.5, line=1)

#dev.off()
